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Abstract 

We introduce a model for describing the dynamics of large numbers 
of interacting cells. The fundamental dynamical variables in the model 
are sub-cellular elements, which interact with each other through phe- 
nomenological intra- and inter-cellular potentials. Advantages of the 
model include i) adaptive cell-shape dynamics, ii) flexible accommo- 
dation of additional intra-cellular biology, and iii) the absence of an 
underlying grid. We present here a detailed description of the model, 
and use successive mean-field approximations to connect it to more 
coarse-grained approaches, such as discrete cell-based algorithms and 
coupled partial differential equations. We also discuss efficient algo- 
rithms for encoding the model, and give an example of a simulation 
of an epithelial sheet. Given the biological flexibility of the model, 
we propose that it can be used effectively for modeling a range of 
multi-cellular processes, such as tumor dynamics and embryogenesis. 
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1 Introduction 



Computational modeling of multi-cellular systems has become an increas- 
ingly useful tool for the interpretation and understanding of experimental 
data in a variety of biological areas. Numerous studies have been performed 
within the modeling community, with application to collective dynamics of 
unicellular organisms, such as myxobacteria ^3] and slime molds ^Tj EHl , 
and to dynamics within multicellular organisms. Studies of the latter in- 
clude avascular tumor growth [SlESl, tumor angiogenesis jlj, embryogenesis 
mini, and cell sorting 

A fundamental issue in modeling cell populations is that of scale. One 
is typically interested in systems composed of tens of thousands to many 
millions of cells, and yet the cell population is often phenotypically heteroge- 
neous. Therein lies the problem of whether or not to include more biological 
realism at the cellular level, with the inevitable computational cost of being 
limited to smaller numbers of cells. The different modeling techniques cur- 
rently employed can be viewed as different compromises to this inescapable 
problem. At the most coarse-grained level one erases cell identity and uses 
continuous cell densities to describe the system. The classic model of this 
type is the Keller-Segel differential equation model of aggregation in social 
amoeba ^^I- Similar differential equation models have been applied to many 
other areas, including, of course, tumor growth jHHHHn]. Cell densities can 
also be modeled using finite element methods jS]. At the next finer scale, 
cells within the population are modeled as discrete objects, yet with little 
or no internal structure [21I2II!- Such models may be constructed either as 
cellular automata on a grid, or else as many-body simulations with no under- 
lying lattice. Recent work has indicated that in the presence of chemotaxis, 
grid effects can lead to strong artifacts Proceeding to smaller length 

scales, cells are endowed with a size (which can change with time during 
the cell cycle), and perhaps anisotropy (e.g. modeled as ellipsoids) to mimic 
cell shape and/or cell polarity j71l221- An ingenious and popular model of 
this type is that due to Graner and Glazier, in which cells are represented 
as clusters of Potts spins on a fine grid TUj. The Potts model approach has 
been applied to a wide range of systems, including cell sorting [Hj, slug for- 
mation in Dictyostelium [TJj, and avascular tumor growth |2H]- At smaller 
scales still, the internal biochemical or biomechanical dynamics of the cell are 
included, e.g. signal transduction for response to chemical signals, and/or 
cytoskeleton dynamics via actin polymerization As the scale of bio- 

logical detail becomes more fine, computational constraints limit the size of 
the cell population. Indeed, models primarily concerned with cytoskeleton 
dynamics typically focus on a single cell. 
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In this paper we introduce a framework for modeling multi-cellular sys- 
tems which is designed to allow simulation of large numbers of cells in three 
dimensions, but also allows for adaptive cell shape dynamics and the ac- 
commodation of successive degrees of intra-cellular biology. This framework 
uses "sub-cellular elements" (defined below) as the fundamental dynamical 
variables, along with overdamped Langevin dynamics [201 12^] for temporal 
development of the system. 

The outline of the paper is as follows. We describe the sub-cellular ele- 
ment model in the next section. In section 3, we use a succession of mean-field 
approximations to connect our model to more coarse-grained descriptions. In 
section 4 we briefly discuss efficient implementation of the model and give a 
simple example of numerical output. We conclude in section 5 with a sum- 
mary of the main results of the paper and a discussion of biological extensions 
of the model. 

2 The sub-cellular element model 

The cell sets the fundamental scale in multi-cellular systems. As such, it is 
natural to base model descriptions of these systems at the cellular scale. One 
of the key properties to incorporate in a cell-based model is dynamical change 
in cell shape (or more generally, cell polarity), which can occur in response 
to local mechanical interactions with neighboring cells, or in response to 
long-ranged chemical signaling pH. Adaptive changes in cell shape/polarity 
allow coherent dynamics of large numbers of cells. For example, several 
mechanisms of large-scale morphological change during gastrulation are due 
to cell intercalation, which is driven by elongation of individual cells along a 
particular axis [21] . From a modeling perspective, cell shape is difficult to pa- 
rameterize. For instance, systematically extending an ellipsoid model of cells 
in three dimensions requires complicated geometrical constructions. Ideally, 
one would like a model in which cell shape emerges from cellular interactions 
- in other words, for cell shape to be adaptive to the local environment. Here, 
we attempt to instantiate this property by sub-dividing each cell into a num- 
ber of "sub-cellular elements." Both the intra- and inter-cellular dynamics 
are written in terms of interactions between these elements. We shall first 
describe this djTiamics by writing the equations of motion for the elements, 
and then discuss how this dynamics can be interfaced with the underlying 
biology. 

For simplicity, consider a system with a constant number N of cells in 
three spatial dimensions, with each cell being composed of M elements. We 
label an individual cell by i G (1, A^) and an element in cell i by ctj G (1, M). 
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For ease of discussion we assume that chemical signahng is absent from the 
system, so that cells respond purely to local biomechanical interactions. In 
this case, the position vector of element is taken to change in time ac- 
cording to three processes: i) a weak stochastic component, which mimics 
the underlying fluctuations in the dynamics of the cellular cytoskeleton, ii) 
an elastic response to intra-cellular biomechanical forces, and iii) an elastic 
response to inter-cellular biomechanical forces. We assume further that the 
elements' motion is over-damped, so that inertial effects can be ignored. The 
equation of motion for the position vector of element takes the form: 

Ya, = Va, - Va, ^ ^ntra(|ya, - Yftl) - Va, ^ ^ Mnter ( | Ya, " Y^jl) ■ (1) 

On the right-hand-side, the noise term t^q,- is a Gaussian distributed random 
variate with zero mean and correlator 

ivZm^t')) = 2vb.,b^^,,fi^-b{t - t') , (2) 

where m and n are vector component labels in the three-dimensional space. 
The second and third terms on the right-hand-side of Eq. (0) represent, re- 
spectively, intra- and inter-cellular interactions between the elements. These 
interactions are completely characterized by the phenomenological poten- 
tials V^ntra and l^nter- At this levcl of description, all relevant biological detail 
must be encoded into these two potentials. The elemental composition of 
cells, along with the inter-elemental potentials, are shown schematically in 
Fig. 1. We have assumed that "two-body" potentials are sufficient to de- 
scribe the dynamics. It may be necessary to use "three-body" potentials to 
capture the essence of more complicated interactions. 

For given biological applications of this modeling framework, one must 
intuit (or better, derive) reasonable forms for Mntra and Mnter- For illustrative 
purposes, consider a population of cells which are weakly adhesive to one an- 
other. Sub-cellular elements both within and between cells will be mutually 
repulsive if their separation is below the equilibrium size of an element. For 
separations larger than this size, the elements will be mutually attractive, 
but with the strength of attraction falling off rapidly with separation. These 
properties can, for example, be conveniently encoded via a generalized form 
of the Morse potential, which is commonly used in physics and chemistry to 
model inter-molecular interactions [25j. The (generalized) Morse potential 
has the explicit form 

V{r) = [/oexp(-r/6) - Voexp{-r/^2) , (3) 

and is illustrated in Fig. 2. It is straightforward to evaluate the position and 
depth of the attractive potential minimum in terms of the four parameters 
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([/q, Vo, ,^1, ^2)- In a simple application of the element model, one can use 
Morse potentials for both l^ntra and Venter, with parameters chosen to ensure 
that the former has stronger inter-elemental adhesion than the latter. This 
condition is necessary, in this simplest version of the model, to maintain the 
mechanical integrity of the cells. 

The introduction and explicit choice of these potentials has so far been 
purely phenomenological, and some discussion of the biological motivation 
for these potentials is necessary. Considering a "typical" tissue cell, such 
as a fibroblast or epithelial cell, the mechanical integrity of the cell is main- 
tained by the internal cytoskeleton P] . This is a complex network of different 
types of interconnected filaments, with actin being the most important fil- 
ament type for cell motility. Dividing the cell into elements corresponds to 
modeling the shape and mechanical integrity of the cell in terms of volume el- 
ements of cytoskeleton. The intra-cellular attraction between elements arises 
from the mechanical rigidity of the cytoskeleton, more specifically, the elas- 
tic forces transmitted through filaments connecting neighboring elements. 
These interactions are local and thus it is necessary that the potential has a 
rapid decay with distance. Elements at opposite sides of a cell mechanically 
interact through elastic forces mediated by elements comprising the interior 
of the cell. The biochemical and biomechanical interactions between cells 
is complex, and arises from a variety of cell-cell (and cell-matrix) contacts, 
such as gap, tight, and anchoring junctions [T]. Still, the interactions are 
local, and, once the cells are linked, one can think of the interaction in terms 
of a short-ranged elastic potential. There is no reason to favor the Morse 
potential at this level of description - there are many reasonable potentials 
that one can write down. Such potentials will, however, be characterized 
by at least four parameters - two energy scales (for short-ranged repulsion, 
and intermediate-range adhesion) and two length scales, which characterize 
the size of an element, and its adhesive range. There are a number of mod- 
els which have been focused on the detailed mechanics of the cytoskeleton, 
and its role in cell motility jHl El UHl An interesting subject of future 
study is the derivation of inter-elemental potentials from coarse-graining the 
underlying cytoskeletal mechanics considered in these more detailed studies. 



3 Connections to coarse-grained models 

In this technical section we sketch the derivation of coarser-grained models 
by applying a succession of mean-field approximations to the element model. 
The essence of this section is summarized in Fig. 3, which shows the funda- 
mental objects/fields characterizing the cell models at different scales. 
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In the first of tliese coarse-graining steps, we replace tlie element model 
by a sub-cellular density model, in which the discrete elements within a given 
cell i are replaced by a smooth average density field Pi(x, t). We stress that a 
separate density field exists for each cell in the system, although the density 
fields are strongly correlated to one another. To proceed, we first recast 
the sub-cellular element model in terms of the probability distribution of 
individual elements. We define the probability distribution of element ctj by 
Pa^(x, t) = (5^(x — YaJ), where the angled brackets denote an average over 
the noise r]. Starting from Eqs. ((H) and we use standard methods piUI^ 
to derive an equation of motion for Pa-, which takes the form 

dtPaA^,t)= Z/V'P„,(x,t) 

+ V- Jd'x' [Vl^ntra(lx-X'l)] ^ P«„^^ (x, t ; x', t) 

+ V- /dV [V\4nter(|x-X'|)]EE^".,/5.(X'^;X''^) '(4) 

jy^i ft- 

where Pa,,i3j is the "two-element" distribution function. The equation of 
motion for this two-element distribution will involve the three-element dis- 
tribution, and so on. The simplest truncation scheme to break the hierarchy 
of equations is the mean field approximation (MFA), in which the statisti- 
cal correlations between elements are discarded. Within this MFA we have 
Pa,,f,^ (x, t; x', t) = Pa, (x, t)Pfs^ (x', t) . 

We now define the sub-cellular density of cell i via Pi(x, t) = J^at PaX'^y t). 
Summing over in Eq. (^, and imposing the MFA, we find a closed equation 
for this sub-cellular density function, which takes the form of an advection- 
diffusion equation: 

dtPii^, t) = uV^Pii^, t) + V ■ pi(x, t) V$,(x, t) , (5) 

where the velocity potential experienced by the density field of cell i is given 
by 

$i(x,t) = J d?x' yintra(|x-x'|)p,(x',t) + j dV Winter ( |x - x' | ) ^ P.' (x, t) • 

(6) 

The MFA used to derive this density equation will typically be good when 
the number of elements used to define the cell is very large. The density 
representation may well be interesting to explore from an analytical stand- 
point. However, it is probably not so useful for numerical implementation. 
For simulation of cells, one must simultaneously integrate N coupled par- 
tial differential equations on a fine three dimensional grid. As we shall see in 
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the next section, the underlying element model, expressed in Eq. (^, can be 
very efficiently encoded for simulation with no need of an underlying grid. 

We take this opportunity to mention that Eq. can be discretized 
in such a way that it resembles a master equation ^T]. In this representa- 
tion, the density of cell i can be interpreted as a probability distribution of 
identical elements, which move from one grid site to the next by activated 
hopping. The hopping rate has the Arrhenius form ~ exp(— A$j/2i/) (where 
A$j is the change in velocity potential, for an element from cell i, between 
the two grid sites of interest, which is actually non-trivial to compute since 
it depends self-consistently on the density pi). The multi-cellular system as 
a whole is then defined on a grid, with each grid site able to accommodate 
(one or more) elements from the different cells. The elements move about 
the lattice (and consequently interact) via activated hopping, with highly 
non-linear hopping rates as indicated above. This representation, although 
not easily implemented, illustrates a qualitative connection between a dis- 
cretized form of the sub-cellular element model (after one level of MFA) 
and the lattice-based Potts model JU] • This discretized form of the element 
model has the flavor of a lattice-gas analog to the Potts model - in the sense 
that a lattice-based element moves over the lattice and yet keeps its original 
parent cell identity, whereas a Potts spin is defined at a lattice site and iden- 
tifies, at a given time, the cell spanning that particular site. Having to hand 
these two distinct models of multi-cellular systems (the sub-cellular element 
model, as expressed in Eq.((TI), and the Potts model) defined at similar scales 
of biological realism, will allow useful cross-validation of these approaches, 
especially when applied to complicated biological systems. 

We can use the density equation to coarse-grain to another scale - 
where now only gross properties (which we refer to loosely as "moments" ) of 
the sub-cellular density field are used to characterize the cell. This coarse- 
graining step is analogous to a multipole expansion in electromagnetism. The 
zeroth moment of cell i is its mass, which is defined by mi{t) = J d^x Pii'x, t). 
Within the present discussion, this quantity is independent of time and cell 
index i since we have assumed that all cells have the same number of elements, 
and that the number of elements does not change with time. Proceeding to 
the ffist moment, we define the position vector of the center of mass of cell 
i via Xj(t) = / d^x xpi(x, t). The equation of motion for this position vector 
is obtained from Eq. ^ and takes the form 



We briefly mention the second moment of cell i, namely its inertia tensor [TH] . 
This is defined via T^'^if) = J d^x {x'^6^^ — x"^x^) pi{-x,t). An equation 
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of motion, similar to Eq. ((7j), can be written down for this tensor. This 
quantity contains crucial information regarding the mechanical polarity of 
the cell. Higher order moments can be defined and contain successively more 
information about the shape and density distribution of the cell. 

Since the equations of motion of these moments are written in terms of 
integrals over the sub-cellular density, they are all strongly inter-dependent. 
To write closed equations again requires some form of truncation. Here we 
use the simplest, which is again a form of MFA, namely: / d^x /(x)pj(x, t) = 
/(xj(t)). This allows us to express the right-hand-side of Eq. (0) in terms of 
Xj(t), and we find the closed equation 

Mt) = -V^V^nter(|Xi(t) -X,(t)|) . (8) 

There are two interesting points to note: i) the intra-cellular potential has 
vanished under this approximation, since we are essentially shrinking the 
cells to points, and ii) the dynamics are now deterministic. Concerning the 
first point, the intra-cellular potential will reappear in this coarse-grained de- 
scription if we include second-order effects - namely, if we derive two coupled 
equations for each cell, describing the time-dependence of the cell's position 
vector and its inertia tensor. Concerning the second point, the effect of the 
noise rja^ has vanished since the first MFA leading to Eq. (0) essentially 
assumes an infinite number of elements, so that the explicit noise terms are 
averaged to zero. The noise from a finite number (M) of elements will be 
non-zero, and have a variance which scales as 1/M. This weak noise (which 
describes the random wandering of the center of mass) can be added to Eq. 
(jHI) a posteriori in order to retain stochasticity in this discrete cell represen- 
tation. One can then write Eq. (jH)) as 

Xi(t) =r/i(t) - V5]Vlntcr(|x.(t) -x,(t)|) , (9) 

where the noise r]i has zero mean, and correlator {f]^{t)ri^{t')) = 2D6ij6"^"'6{t— 
t'), where D = u/M. This stochastic model, which tracks the positions of 
the cells, is precisely that studied by Newman and Grima As shown in 
that work, a further MFA applied to Eq. leads to a closed equation for 
the density of cells, which is defined via n(x, 

t) = Ei('^^(x-x,(t))). We omit 
the details here and simply give the final result: 

dtn{x, t) = L)V^n(x, t) + V ■ n(x, t)V^(x, t) , (10) 

where the coarse-grained velocity potential \I' for the cell density has the 
form 

^(X,t) = I Sx' Vi,ter(|x-X'|)n(x',t) . (11) 
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A rigorous derivation of this last step has been given by Stevens in the 
context of chemotaxis j^ZI, and uses the hmit of infinite cell number, with 
an appropriately scaled chemotactic coupling. 

Finally, after three levels of coarse-graining, we have arrived at a partial 
differential equation for the cell density, as given in Eqs. (fTnjl and (fTT|l . 
As mentioned in the Introduction, this level of description has been widely 
used to describe the large-scale dynamics of cell populations. However, as 
should be clear from this analysis, a great deal of statistical information 
and smaller-scale biomechanics must be discarded at this scale. It would 
be very interesting to rederive the density equations from a more careful 
analysis. Some details of the intra-cellular potentials (especially regarding 
cell polarity) can be captured in this largest-scale description through i) 
calculating renormalized parameters, such as the diffusion coefficient D, in 
the density description (|Tnjl . and ii) deriving the companion equation for a 
"cell polarity field" from the discrete cell equations for the inertia tensor. 

4 Efficient algorithms and model output 

We now return to the sub-cellular element model as described in section 2. 
The numerical implementation of this model turns out to be fairly straight- 
forward and efficient. Since the fundamental dynamical variables are position 
vectors, we have no need for an underlying grid, and simply need to track 
the values of the (M x N) vectors {ya,} which completely describe the state 
of the system at any given time. 

It is worth mentioning that some care must be taken in constructing the 
algorithm so as to avoid a CPU cost which scales as (MN)"^. This would arise 
from attempting to interact every element with every other in order to update 
the system. Clearly this is not necessary since the potentials are short-ranged. 
As such, the algorithm only needs to interact a given element with those 
elements which are close enough to have a non-negligible interaction. So long 
as we can efficiently identify these nearby elements, our algorithm will have 
a CPU cost which scales as MN. This will allow the simulation of large 
numbers of cells with moderate to large numbers of elements per cell. There 
are a number of ways to identify nearby elements. The methods to achieve 
this have been developed over the years in molecular dynamics simulations 
1^ . Examples are neighbor tables and the more sophisticated binary 
search trees and octrees. We have employed a method based on "sectors." 
The three dimensional system is broken up into a grid of sectors, the size of a 
sector chosen to be about twice the range of inter-elemental interactions. The 
dynamics of the elements are completely oblivious to the sectors. The sectors 
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simply allow one to construct a look-up table wherein for each sector there 
is a list of the identity of those elements in that sector. When one calculates 
the interactions of a given element, one computes the interactions between 
that element and all those in its own sector and the nearest-neighbor sectors. 
Temporal development is performed using an explicit Euler discretization, 
with time increment 5t. 

In Fig. 4 we show an example of the output of this algorithm. In this 
example we have simulated 128 cells, with each cell constructed from 20 
sub-cellular elements. Both the intra- and inter-cellular potentials are cho- 
sen to be generalized Morse potentials, with parameter sets (f/o, Vo,.Ci,^2) = 
(0.25, 0.1, 0.12, 0.36) and (0.25, 0.05, 0.12, 0.24) respectively. Other parameter 
values are v = 0.001 and St = 0.1. For ease of presentation we have simu- 
lated the cells in a quasi-two-dimensional geometry, with hard-wall boundary 
conditions at 2; = and z = 0.5, which could represent an epithelial layer 
bounded by a basement membrane. The system is shown after about 3000 
iterations, starting from a random distribution of cell positions (with the 
initial positions of the elements for each cell randomly distributed in a small 
region about these cell positions). This simulation requires about 3 minutes 
on a 2GHz PC. Extrapolating to larger systems, we see that thousands of 
iterations for a system of 10, 000 cells (each with 20 elements) requires a few 
hours of CPU time on a PC. More sophisticated optimization of the algo- 
rithm will allow further improvements in efficiency. Note that as a result of 
the elements attempting to form inter-cellular adhesive bonds with elements 
of nearby cells, the cells have adapted their cell shapes to their local biome- 
chanical environment. A more systematic numerical study comparing the 
different coarse-grained descriptions is currently in progress (26j . 

5 Summary and outlook 

In this paper we have introduced a model of interacting multi-cellular sys- 
tems, in which the fundamental objects are not cells, but "sub-cellular ele- 
ments" - by which we mean, in the simplest sense, small volume elements 
of intra-cellular cytoskeleton. The dynamics of the elements is described by 
Langevin equations, as given in Eq. (0). The three dynamical contribu- 
tions to a given element's motion are i) a weak stochastic component, ii) 
local biomechanical interactions with other elements within the same cell, 
described by a phenomenological potential V^ntra, and iii) local biomechani- 
cal interactions with elements in nearby cells, these described by a potential 
Venter (sce Fig. 1). The success of the model in a given biological application 
depends, in large part, on well-chosen forms of these potentials. The generic 
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form of these potentials is reasonably well captured by the generalized Morse 
potential (Fig. 2). We have outlined, in section 3, a series of coarse-graining 
procedures (summarized in Fig. 3), whereby the sub-cellular element model 
can be linked to models at larger scales, such as discrete cell models and dif- 
ferential equation models describing the macroscopic cell density. In section 
4 we indicated how an efficient algorithm may be constructed to integrate 
Eq. (PJ) forward in time for large numbers of cells/elements, and gave an 
example from a simulation of a sheet of cells (Fig. 4). 

We end the paper with a discussion of some of the many possible exten- 
sions of the model, whereby biological detail can be added without distortion 
of the underlying element framework. 



For simplicity, in this introductory paper we have presented the element 
model for a population of N identical cells. Phenotypic heterogeneity at 
the cell level can be described by attaching cell labels to the noise strength 
(z/ z/j), intra-cellular potentials (Vintra Vi), and inter-cellular potentials 



5.2 Sub-cellular element types 

For many applications it will not be sufficient to construct a cell from M 
identical elements. For instance, it might be necessary to distinguish el- 
ements in the interior of the cell ("cytoplasmic elements") from elements 
on the surface of the cell ("membrane elements"). Different element types 
can easily be instantiated by defining the appropriate classes of intra- and 
inter-cellular potentials. In the example just given, the inter-cellular interac- 
tions would primarily be described by an inter-cellular potential connecting 
membrane elements from neighboring cells, while the cytoplasmic elements 
would interact only with elements within the same cell. In addition, elements 
can be endowed with internal variables registering environmental variables 
such as pH or nutrient level, and communicate these variables to neighboring 
elements to trigger appropriate cell response. 

5.3 Extra- cellular element types 

Cells not only adhere to each other, but also interact biomechanically with 
a range of non-cellular environmental structures, such as the extra- cellular 
matrix or various gel-like media [T]. In some cases these structures are pro- 
duced by the cells themselves. Within the spirit of the element model, these 



5.1 Cell types 
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extra-cellular structures can also be constructed from elements (i.e. elasti- 
cally coupled degrees of freedom on the same length scale as the sub-cellular 
elements) with appropriately defined interaction potentials to describe the 
cell/non-cell interactions. 

5.4 Chemotaxis 

We have been exclusively concerned with short-ranged biomechanical inter- 
actions in this paper. Equally important for many applications (e.g. embryo- 
genesis, wound healing) are long-ranged chemical interactions. The simplest 
way to introduce such interactions in the element model is to chemically cou- 
ple the centers of mass of cells which are signaling to one another. Then, since 
the source and sink of the signal are points, it is straightforward to imple- 
ment the Green function methods developed in some detail in Newman and 
Grima jSD]- These methods encode the diffusion of chemical signals through 
the diffusion equation Green function (which allows one to avoid introduc- 
ing a fine grid for explicit integration of the chemical diffusion equations). 
More sophisticated treatments would be based on cell polarity induced by the 
chemical signal. For instance, one could use chemically responsive element 
types within the cell, so that the cell as a whole responds to a chemical signal 
via a sub-cellular element response, followed by a cell-level response mediated 
by elastic interactions of the responsive element to the non-responsive ones. 

5.5 Cell cycle 

For many applications, and in particular that of tumor growth, cells in the 
population are undergoing numerous cell divisions over the time scales of 
interest. Cell growth can be accommodated in the modeling framework by 
allowing new elements to be spawned within the cell (based on, for instance, 
an internal monitoring of nutrient levels as discussed above in subsection 
5.2). Mitosis is more complicated. It can be "forced" upon a cell using some 
form of threshold conditions, under which the elements segregate into two 
daughter cells. However, more elegant mechanisms can no doubt be devised. 

In this final section we have tried to give a flavor of possible biological 
extensions of the model, and the ease with which they can be implemented 
within the element framework. This flexibility must ideally be tempered by 
minimal incremental model-building in order to keep models simple enough 
to understand from both quantitative and biological perspectives. With this 
caveat in mind, we hope that the sub-cellular element model will prove useful 
in the computational study of a wide range of multi-cellular systems. 



12 



The author would hke to thank James Glazier, Ramon Grima, John Nagy, 
Kevin Schmidt, Erick Smith, and Cornelius Weijer for interesting discussions. 
The author gratefully acknowledges partial support from the NSF (lOB- 
0540680) and NSF/NIH (DMS/NIGMS-0342388). 



13 



References 



[I] B. Alberts, A. Johnson, J. Lewis, M. Raff, K. Roberts, and P. Walter, 
Molecular biology of the cell, 4th edition. Garland Publishing, New York, 
2004. 

[2] T. Bretschncidcr, B. Vasiev, and C. J. Wcijcr, A model for Dictyostelium 
slug movement, J. Theor. Biol. 199 (1999), 125-136. 

[3] H. M. Byrne, Modeling avascular tumor growth, in Cancer Modelling and 
Simulation, L. Preziosi, ed.. Chapman and Hall/CRC, Boca Raton, FL, 
2003. 

[4] M. A. J. Chaplain, Mathematical modeling of angiogenesis, J. Neuroon- 
col. 50 (2000), 37-51. 

[5] L. A. Davidson, M. A. R. Koehl, R. Keller, and G. F. Oster, How do sea 
urchins invaginate - using biomechanics to distinguish between mecha- 
nisms of primary invagination. Development 121 (1995), 2005-2018. 

[6] P. A. Dimilla, K. Barbee, and D. A. Lauffenburger, Mathematical model 
for the effect of adhesive mechanics on cell-migration speed, Biophys. J. 
60 (1991), 15-37. 

[7] D. Drasdo, R. Krcc, and J. S. McCaskill, Monte Carlo approach to tissue- 
cell populations, Phys. Rev. E 52 (1995), 6635-6657. 

[8] J. A. Glazier, and F. Graner, Simulation of the differential adhesion 
driven rearrangement of biological cells, Phys. Rev. E 47 (1993), 2128- 
2154. 

[9] M. E. Gracheva and H. G. Othmer, A continuum model of motility in 
amoeboid cells. Bull. Math. Biol. 66 (2004), 167-193. 

[10] F. Graner, and J. A. Glazier, Simulation of biological cell-sorting using 
a two-dimensional extended Potts model, Phys. Rev. Lett. 69 (1993), 
2013-2016. 

[II] R. Grima, and T. J. Newman, Accurate discretization of advection- 
diffusion equations, Phys. Rev. E 70 (2004), #036703. 

[12] R. Grima, S. Sa-nguansin, and T. J. Newman, in preparation (2005). 

[13] J. M. Haile, Molecular dynamics simulation: elementary methods, Wiley 
Interscience, 1997. 



14 



[14] O. A. Igoshin, A. Mogilner, R. D. Welch, D. Kaiser, and G. Oster, Pat- 
tern formation and traveling waves in myxobacteria: theory and model- 
ing, Proc. Natl. Acad. Sci. (USA) 98 (2001), 14013-14018. 

[15] E. F. Keller, and L. A. Segel, Initiation of slime mold aggregation viewed 
as an instability, J. Theor. Biol. 26 (1970), 399-415. 

[16] L. D. Landau and E. M. Lifshitz, Mechanics, 3rd ed., Butterworth- 
Heinemann, Oxford, 2000. 

[17] A. F. M. Maree, and P. Hogeweg, How amoeboids self-organize into 
a fruiting body: multicellular coordination in Dictyostelium discoideum, 
Proc. Natl. Acad. Sci. USA 98 (2001), 3879-3883. 

[18] A. Mogilner and G. Oster, Cell motility driven by actin polymerization, 
Biophys. J. 71 (1996), 3030-3045. 

[19] J. Nagy, The ecology and evolutionary biology of cancer: a review of 
mathematical models of necrosis and tumor cell diversity. Math. Biosci. 
Eng. 2 (2005), 381-418. 

[20] T. J. Newman, and R. Grima, Many-body theory of chemotactic cell-cell 
interactions, Phys. Rev. E 70 (2004) #051916. 

[21] K. J. Painter, P. K. Maini, and H. G. Othmer, A chemotactic model 
for the advance and retreat of the primitive streak in avian development. 
Bull. Math. Biol. 62 (2000), 501-525. 

[22] E. Palsson and H. G. Othmer, A model for individual and collective cell 
movement in Dictyostelium discoideum, Proc. Natl. Acad. Sci. (USA) 97 
(2000), 10448-10453. 

[23] D. C. Rapaport, The art of molecular dynamics simulation, Cambridge 
University Press, 2004. 

[24] B. Rubenstein, K. Jacobson, and A. Mogilner, Multiscale two- 
dimensional modeling of a motile simple-shaped cell, Multiscale Model. 
Sim. 3 (2005), 413-439. 

[25] L. I. Schiff, Quantum mechanics, 3rded., McGraw-Hill, Singapore, 1968. 

[26] E. Smith, T. J. Newman, and C. J. Weijer, in preparation (2005). 

[27] A. Stevens, The derivation of chcmotaxis equations as a limit dynam- 
ics of moderately interacting stochastic many-particle systems, SIAM J. 
Appl. Math. 61 (2000), 183-212. 



15 



[28] E. L. Stott, N. F. Britton, J. A. Glazier, and M. Zajac, Stochastic sim- 
ulation of benign avascular tumor growth using the Potts model, Math. 
Comp. Model. 30 (1999), 183-198. 

[29] N. G. van Kampen, Stochastic processes in physics and chemistry. North 
Holland, Amsterdam, 1992. 

[30] B. Vasiev and C. J. Weijer, Modeling of Dictyostelium discoideum slug 
migration, J. Theor. Biol. 223 (2003), 347-359. 

[31] L. Wolpert, R. Beddington, T. Jessell, P. Lawrence, E. Meyerowitz, and 
J. Smith, Principles of development, 2nd ed., Oxford University Press, 
2002. 



16 




Figure 1: Schematic diagram showing two cells, i and j, and a subset of the 
intra- and inter- cellular interactions between their elements. The elements 
of cell i are represented by open circles, and those of cell j by filled circles. 
The intra- and inter-cellular interactions are represented by solid and dashed 
lines respectively. 



17 



0.16 




-0.04 I ^ ^ ^ ^ 1 

0.2 0.4 0.6 0.8 1 



Figure 2: The Morse potential for parameter values Uq = 0.25. Vq = 0.1, 
^1 = 0.12, and ^2 — 0.36, as used for the intra-cellular potential in section 4. 
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Figure 3: Relationships between model descriptions of multi-cellular systems 
at different scales. The four levels (from finer to coarser scales) are described 
explicitly by Eqs. ((H), (jSj), ©, and (fTIH) respectively. 
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Figure 4: An example of numerical output from the model using Morse 
potentials. Parameter values are given in the main text. Shown here is 

a two-dimensional projection of data from a model of an epithelial sheet. 
Each of the 128 cells is composed of 20 elements (open circles), and element 
motion is constrained in the third dimension by hard-wall boundaries. The 
filled circles indicate the center of mass of each cell. 
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